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ABSTRACT 


The  objective  of  this  investigation  is  to  study  the  damage  mechanics  of  composite 
structures  using  a  micromechanical  approach  for  determining  strength  and  stiffness 
degradation  of  the  composite  structures  as  damages,  such  as  matrix  cracking  and  fiber 
breakage,  progress.  The  micromechanical  cell  method  provides  for  analysis  of  stress  at 
the  fiber  and  matrix  level  while  providing  smeared  composite  properties  for  global 
structural  analysis.  As  a  result,  the  damage  and  failure  criteria  are  expressed  in  terms  of 
the  fiber  and  matrix  stress  level  of  the  composite  structure.  A  correlation  for  stiffness 
reduction  due  to  transverse  cracking  of  a  ceramic  matrix  composite  under  tensile  loading 
is  implemented  in  a  three-dimensional  finite  element  model.  Next,  thermal  residual 
stresses  from  fabrication  of  the  ceramic  matrix  composite  are  incorporated  into  the 
analysis.  Finally,  the  finite  element  method  is  applied  to  a  polymer  matrix  composite 
laminate  with  a  center  hole  in  order  to  study  the  progression  of  damage  and  final  failure 
during  tensile  loading.  The  comparisons  between  the  present  predictions  and  the 
experimental  results  for  the  previous  examples  are  very  good. 
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L  INTRODUCTION 


A.  BACKGROUND 

Modem  structural  designs  increasingly  incorporate  man-made  composite  materials 
in  applications  that  require  components  with  special  material  properties  which  are 
unavailable  from  conventional  metals  or  alloys.  From  the  structural  mechanics  viewpoint, 
composites  are  typically  used  to  improve  either  the  stiffness-to-weight  ratio  or  the 
strength-to-weight  ratio  of  a  structural  member.  Naval  shipboard  applications  of 
composite  materials  are  typically  dependent  on  these  weight  savings  when  a  composite 
material  is  used  to  fulfill  a  specified  stiffness  or  strength  requirement.  Current  naval 
composite  applications  include  main  propulsion  shafting,  decks,  superstructures,  small 
equipment  foundations  and  graphite-reinforced-poly  (GRP)  auxiliary  system  piping. 
Although  the  composite  material  itself  may  actually  cost  more  than  its  conventional 
metallic  predecessor,  incorporation  of  lightweight  composite  materials  usually  reduces 
total  ship  cost  due  to  synergistic  savings  in  the  entire  ship  structure  because  of  the 
reduced  total  ship  weight. 

Composite  properties  such  as  strength  and  stiffness  are  dependent  upon  the  volume 
fraction  of  the  fiber  and  the  individual  properties  of  the  constituent  fiber  and  matrix 
materials.  The  variation  of  the  fiber  volume  fraction  for  a  given  composite  allows  the 
designer  greater  flexibility  when  incorporating  composite  materials  into  a  structure.  With 
this  flexibility  comes  added  complexity  due  to  the  fact  that  the  optimum  volume  fraction 
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for  a  given  stiffness  may  not  be  optimum  for  composite  strength.  And,  composite 
material  properties  may  vary  due  to  damage  accumulation  such  as  matrix  transverse 
cracking  during  the  loading  history  of  the  composite  member.  This  current  research 
investigates  the  analysis  of  composite  material  properties,  such  as  stiffness  and  strength, 
and  their  degradation  using  a  three  dimensional  micromechanical  model. 

B.  DETERMINING  COMPOSITE  MATERIAL  PROPERTIES 

The  study  of  composite  materials  and  structures  can  be  undertaken  from  two 
different  approaches;  micromechanical  and  macromechanical.  In  the  micromechanical 
approach,  the  properties  of  the  constituent  fiber  and  matrix  materials  and  their  interaction 
through  stress-strain  constitutive  relations  are  analyzed  in  order  to  predict  the  overall 
behavior  of  the  composite  structural  member.  In  the  macromechanical  approach,  the 
properties  of  the  constituent  fiber  and  matrix  materials  are  averaged  or  smeared  to 
produce  a  set  of  pseudo-homogeneous  properties  for  the  composite  structural  member. 
The  macromechanical  approach  has  the  advantage  of  requiring  less  detailed  modeling  in 
that  individual  fiber  and  matrix  properties  are  only  used  initially  in  determining  the 
smeared  composite  properties.  Finite  element  models  using  the  macromechanical 
approach  are  somewhat  less  complex  and  typically  require  less  computational  time  than 
those  models  v  hich  are  based  on  the  micromechanical  approach.  However,  the 
micromechanical  approach  provides  more  useful  information  such  as  the  fiber  and  matrix 
stress  and  strain  which  are  typically  used  for  failure  criteria  or  strength  degradation 
computations.  Thus,  the  more  complex  micromechanical  approach  is  more  useful  when 
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analyzing  damage  mechanics  and  failure  of  composite  material  structures. 

Under  normal,  non-damaging  or  elastic  loading  conditions,  the  stiffness  properties 
of  aligned,  continuous,  fiber  composites  can  be  predicted  by  the  simple  rule  of  mixtures 
(ROM).  The  ROM  model  uses  the  strength  of  materials  approach  and  is  based  upon 
assumptions  that  the  composite  is  loaded  at  low  strain  levels  in  the  elastic  region  where 
no  damage  occurs  in  the  fibers  or  matrix,  that  the  fibers  have  uniform  properties  and  are 
aligned  parallel  throughout  the  composite,  that  the  matrix  and  fibers  have  a  no-slip 
perfectly  bonded  interface,  and  that,  for  the  longitudinally  applied  load  case,  there  is  equal 
strain  in  the  fiber  and  matrix  [Refl].  For  composites  which  approximate  these 
assumptions  and  which  are  loaded  in  the  longitudinal  direction  along  the  fiber,  ROM 
provides  reasonable  values  for  composite  stiffness.  However,  for  transverse  loading  with 
respect  to  the  fiber  axis,  ROM  does  not  accurately  predict  the  composite  stiffness. 

Micromechanical  models  based  on  the  theory  of  elasticity  provide  better  overall 
results  for  aligned,  parallel  fiber  composites  [Ref  2].  Noteworthy  among  the 
micromechanical  models  are  the  relations  developed  by  C.C.  Chamis  [Ref  3]:  Chamis 
assumed  that  the  matrix  is  isotropic  and  that  the  fibers  are  orthotropic.  A  composite  with 
such  constituents  would  have  one  plane  of  symmetry  and  therefore  would  be  transversely 
isotropic.  With  the  condition  of  transverse  isotropy,  Chamis  reduced  the  problem  from  one 
with  three  dimensions  to  one  with  two  dimensions.  This  method  resulted  in  five 
independent  properties  in  order  to  define  a  stress-strain  constitutive  relation  for  the 
smeared  composite.  Chamis  and  Sendeckyj  [Ref  4]  provide  an  excellent  survey  of  the 
early  micromechanics  approach  to  predicting  composite  stiffness.  Jones  [Ref  5] 
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succinctly  presents  methods  for  bounding  stiffness  properties  of  composites  using 
variational  energy  approaches  with  classical  elasticity  theory  and  also  describes  the 
contiguity  approach  used  to  develop  the  famous  HaJpin-Tsai  relations.  These  previous 
methods  provide  generally  reasonable  values  for  composite  stiffness.  But,  these  early 
micromechanical  methods  are  often  not  sufficiently  accurate  or  readily  applied  to 
structural  computations  involving  composites.  And,  these  methods  are  based  on  two- 
dimensional,  linear  elastic  analysis  while  the  real  composite  member  is  three  dimensional, 
and  typically  behaves  in  nonlinear  manner  over  its  useful  loading  range  [Ref  6  ]. 

In  recent  years,  there  has  been  significant  study  of  three-dimensional 
micromechanical  models.  Dvorak  and  Bahei-EI-Din  [Ref  7]  formulated  an  axisymmetric 
model  for  an  elastic-plastic  constitutive  relation  based  on  cylindrical  fibers  with 
surrounding  matrix  material.  The  axisymmetric  assumption  was  considered  physically 
valid  and  reduced  the  complexity  of  the  three  dimensional  elasticity  analysis  for  the 
micromechanical  modelling  of  elastic-plastic  deformation.  Aboudi  [Ref  8]  introduced  the 
micromechanical  method  of  cells  based  on  the  assumption  that  the  composite  material 
consisted  of  two  repetitive  phases  of  fiber  and  matrix  materials.  By  assuming  the 
repetitive  or  periodic  nature  of  the  fiber  and  matrix  array,  Aboudi  was  able  to  simplify 
his  model  by  using  a  representative  cell  which  consisted  of  four  subcells.  Although  a 
view  of  his  model  appears  to  desc  ibe  a  square  cell  with  one  subcell  of  fiber  and  three 
subcells  of  matrix,  the  geometry  of  the  fiber  and  matrix  is  not  restricted  due  to  calculation 
of  the  interface  conditions  on  an  average  bases.  Aboudi  introduced  the  method  of  cells 
starting  with  simple  unidirectional  fiber  composites  and  then  extended  his  method  to 
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discontinuous  short  fiber  composites.  Aboudi's  micromechanical  model  is  mathematically 
complex  and  computational  intensive.  Others,  such  as  Pecknold  [Ref  9]  noted  that 
Aboudi's  model  forms  the  basis  for  a  finite  element  model.  Pecknold  conducted  an 
investigation  of  a  simplified  unit  cell  model.  Kwon  used  a  micromechanics  model  [Ref 
10,  II,  12]  by  focusing  on  the  fiber  and  matrix  stresses  at  the  micromechanical  level. 
Kwon  later  refined  his  original  cell  method  [Ref  6].  Kwon's  micromechanics  models  are 
especially  applicable  for  the  investigation  of  composite  damage  because  it  considers  both 
the  fiber  and  matrix  stress  at  the  micromechanical  level  and  thus  allows  specific  fiber  and 
matrix  yield  and/or  failure  criteria  to  be  applied.  The  refined  Kwon  model  [Ref  6] 
formed  the  basis  for  this  current  work. 

C.  THE  MICROMECHANICAL  CELL  MODEL 

Similar  to  the  Aboudi  cell  model,  the  Kwon  model  considers  the  composite  as  a  unit 
cell  composed  of  four  subcells:  one  fiber  subcell  and  three  matrix  subcells.  The  unit  cell 
is  represented  by  a  three  dimensional  solid,  the  rectangular  parallelpiped.  Based  on 
symmetry,  only  one-quarter  of  the  total  cell  may  be  modeled  as  shown  in  Figure  1.1. 
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Note  that  the  size  of  each  cell  is  dependent  on  the  fiber  volume  fraction.  Kwon 
then  expresses  the  composite  stresses  and  strains  as  a  function  of  the  subcell  stresses,  the 


subcell  strains,  and  the  volume  fraction  as  shown  in  Equations  (1.1)  and  (1.2); 

5,,  -  V^-  (l  (l  *(l  u  .  1,3  (1.1) 

?  =  (l  (l  *  (l  -f^)  8*  i,j  .  1,3  (1.2) 

where  o,j  and  e,j  are  composite  stresses  and  strains,  o“j  and  e“j  are  subcell  (a  =  a,  b,  c, 
or  d)  stresses  and  strains,  and  Vf  is  the  fiber  volume  fraction. 

Note  that  subcell  'a'  represents  the  fiber  and  subcells  'b',  'c',  and 'd'  represent  the 

matrix.  Thus,  the  composite  smeared  stress  or  strain  values  are  equal  to  the  volume 

average  of  the  subcell  stresses  and  strains.  The  condition  of  stress  continuity  is  satisfied 

at  the  subcell  interfaces  as  expressed  by  Equation  (1.3). 
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Each  subcell  may  have  different  strains,  but  the  following  strain  compatibilities  are 
assumed.  The  longitudinal  strains  of  the  subcells  are  equal  and  the  sums  of  the  transverse 
strains  of  the  subcell  in  either  transverse  direction  are  also  equal.  The  shearing  strains 
must  also  satisfy  similar  conditions.  Expressed  mathematically,  the  strain  conditions  are 


shown  in  Equation  (1.4). 
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The  constitutive  relation  for  each  subceil  is  expressed  by  the  generalized  Hooke's 
law  of  Equation  (l  .S). 

=E“k,eu  i,j,k,l  =1,2,3  and  a  =a,b,c,d  (15) 

Solving  Equations  (1.1)  to  (1.5)  together  yields  the  expression  for  smeared 
composite  material  properties  in  terms  of  fiber  and  matrix  material  properties  and  their 
volume  fractions.  In  addition,  the  equations  provide  the  following  sequential  relations; 
global  displacements  — composite  strains  — >  fiber  and  matrix  strains  — >  fiber  and 
matrix  stresses  — ♦  composite  stresses. 

The  global  displacements  are  obtained  from  the  finite  element  analyses  of  fibrous 
composite  structures.  For  damage  progression  studies,  the  calculation  of  stress  and  strain 
at  the  fiber  and  matrix  level  allows  failure  criteria  to  be  applied  at  that  micro-level.  Thus, 
numerical  modeling  using  this  method  can  provide  insight  into  the  micro-level  failure 
process  for  composite  materials. 

For  this  current  work,  the  Kwon  micromechanical  method  of  cells  was  implemented 
within  the  finite  element  program  in  order  to  perform  several  function ::  First,  the  new 
micromechanical  model  was  used  to  determine  the  smeared  composite  properties  for  use 
in  computing  the  basic  finite  element  stiffness  matrix.  Secondly,  the  model  was 
implemented  in  the  finite  element  post-processing  routines  to  convert  the  finite  element 
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displacement  results  into  stresses  and  strains  at  the  fiber  and  matrix  micro-level.  The 
fiber  and  matrix  stress  levels  were  used  for  either  stiffness  reduction  or  the  application 
of  failure  criteria. 
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IL  FINITE  ELEMENT  MODEL  DEVELOPMENT 


Finite  element  numerical  analysis  was  used  to  solve  the  three  dimensional  elasticity 
problem  in  order  to  determine  stress  and  strain  of  the  fiber,  matrix  or  the  overall 
composite  specimen.  The  finite  element  method  transforms  the  partial  differential 
equations  to  a  system  of  algebraic  equations.  The  primary  solution  was  obtained  in  terms 
of  displacement.  A  solution  post-processor  subroutine  converted  the  displacement  solution 
into  composite  strain.  The  Kwon  micromechanical  method  of  cells  provided  the  algorithm 
to  solve  for  local  cell  strains,  local  celt  stresses  and  composite  stress.  Stress  values 
determined  at  the  micromechanical  cell  level  provided  the  dependent  variable  for  stiffness 
reduction  or  failure  relations.  Details  of  the  finite  element  method  derivation  are  given 
below. 

A.  ANALYTICAL  DERIVATION 

The  derivation  for  the  force  equilibrium  equations  of  a  three  dimensional  solid  body 
which  is  experiencing  negligible  body  forces  is  found  in  any  advanced  solid  mechanics 
or  elasticity  textbook  [Ref  13,  14,  15,  16],  These  equations  are  shown  in  Equation  (2.1). 
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These  three  equilibrium  equations  are  written  in  terms  on  nine  stress  variables.  However, 
only  six  of  the  stress  variables  are  independent  due  to  the  requirement  for  moment 
equilibrium.  Application  of  moment  equilibrium  conditions  on  the  unit  solid  element 
yields  Equation  (2.2). 


The  unit 
establish 


=v 

=x_ 


(2.2) 


solid  element  is  shown  in  Figure  2.1  below  with  stress  terms  as  indicated  to 
a  reference  for  the  sign  convention  and  notation. 


Figure  2. 1  Unit  Solid  Cell  With  Stress  Components 


II 


Let  'u',  'v',  and  'w'  represent  displacement  in  the  'x',  'y'.  'z'  directions 

respectively.  Then,  the  relationship  of  strain  to  displacement,  assuming  small 
displacements,  can  be  written  as  shown  in  Equation  (2.3). 
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X 


Y 

»  xy 


5u  _  ^  ^ 

5x  ’  ^  dy  '  *  5z 

=  V  V 

(ay  axj’  dz}'  [dz  axj 


(2.3) 


The  strain-displacement  relations  introduce  six  strain  terms  as  unknown  variables.  Also, 
at  this  point,  'u',  'v',  and  'w'  are  unknown.  These  equations  can  be  written  more 
conveniently  in  matrix  form  as  in  Equation  (2.4). 
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(2.4) 


Next,  consider  the  constitutive  relationship  between  stress  and  strain  of  Equation  (2.S). 


{a)=[D]{8} 


(2.5) 
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In  Equation  (2.5),  {a}  is  a  '6  x  1’  column  vector  of  stress,  {&}  is  a  '6  x  1'  column  vector 
of  strain,  and  [D]  is  the  *6  x  6'  matrix  of  material  properties.  Equations  2.1,  2.3,  and  2.5 
combine  to  give  IS  equations  with  IS  unknown  variables.  Unfortunately,  the  equations 
include  several  partial  differential  equations. 

The  derivation  of  the  finite  element  equations  from  the  three  dimensional  elasticity 
equation  above  is  based  on  course  notes  and  commonly  available  current  textbooks  [Ref 
17,  18,  19,  20].  The  method  of  weighted  residuals  as  modified  by  Galerkin  is  employed 
to  develop  the  displacement  based  finite  element  equations.  Eight  noded  isoparametric 
rectangular  paralleipiped  elements  are  used  for  the  formulation.  Linear  shape  functions 
are  used  for  analytical  and  computational  simplicity.  Detailed  derivation  according  to 
these  precepts  is  continued  below. 

The  first  step  in  converting  the  partial  differential  equations  into  algebraic  equations 
is  to  apply  the  method  of  weighted  residuals  to  the  three  dimensional  stress  equilibrium 
Equation  (2.1).  To  this  end,  each  of  the  three  equations  of  Equation  (2.1)  are  multiplied 
by  a  weight  function,  (i=l,  2,  3),  which  is  continuous  over  the  physical  domain  of  the 
problem.  Then,  the  three  new  product  equations  are  integrated  over  the  entire  problem 
domain.  The  goal  is  to  chose  weight  functions  W,  which  are  orthogonal  to  the  initial 
residuals  of  the  equilibrium  equations  such  that  the  integral  becomes  equal  to  zero.  If 
'V  is  the  domain  volume  of  the  problem,  the  weighted  residual  equations  are  now  shown 
in  Equation  (2.6). 
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At  this  point  it  should  be  noted  that  boundary  conditions  must  be  specified  for  Equation 
(2.1).  The  boundary  conditions  may  be  specified  as  (1)  essential  or  geometric  boundary 
conditions  where  some  surface  displacements  are  specified,  (2)  natural  or  stress  boundary 
conditions  where  surface  tractions  are  specified,  or  (3)  a  combination  of  these  types  of 
boundary  conditions.  Before  applying  boundary  conditions,  further  manipulation  of  the 
weighted  residuals  is  required  as  shown  by  Equation  (2.7);  Integration  by  parts, 
commonly  referred  to  as  Gauss'  theorem  or  the  divergence  theorem  when  applied  in  three 
dimensions,  will  be  performed  in  order  to  weaken  the  continuity  requirements  on  the 
approximate  solution. 


'•  ■  ///  [-('’•5  •  *  // 

''  ■  ///  [■('»?  *  *  //  *  V."'")®  ^ 


In  Equation  (2.7),  'S'  is  the  domain  surface  boundary,  'V  is  the  domain  volume,  and  'n„', 
'ny',  and  'n^'  are  outward  unit  surface  normal  directions  cosines  in  the  'x',  'y',  and  'z' 
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directions,  respectively.  Now,  define  the  boundary  stress  conditions  by  surface  tractions 
terms  as  shown  in  Equation  (2.8); 

<b  =  o  n  +T  n  +T  n 

▼x  X  X  xy  y  xx  x 

(2.8) 

Equation  (2.7)  can  be  written  in  matrix  form  with  Equation  (2.8)  incorporated; 


(W,  SfV.  dfV. 

o, +  — -  +  T„ — - 

*  ac  ^  dj  * 

ttt  J  dW^ 

dW.  cW. 

T_ — -  +  r„ — -  +  a. — - 
“  ^  ^  ^  dL 


dV-fjUf,  dS 


Equation  (2.9)  can  be  further  modified  by  separating  the  column  vector  inside  the  volume 
integral  into  a  product  of  a  '3  x  6’  matrix  and  a  '6  x  1’  vector.  This  step,  shown  in 
equation  (2.10),  isolates  the  weighting  function  derivatives  in  the  matrix  and  the  stress 
terms  in  the  vector. 


d^.  d9f,  dW. 

a - -  +  r  - -  +  r  — - 

'A  ^  df  *  4 

dW.  dW. 

T  - -  +  O  - -  +  T  - - 

^  A  ^  ds  ^  at. 

dfV,  dfVt  Af, 

r  — i  ♦  r  — -  *  a, — - 
*  A  ^  df  *  A 


0  0 


0  0 


m.  aw. 

►  — I  — ?  0 
A  A 

V,  aw,  aw, 
J.  0  — -  — ^ 
t  A'  * 


(2.10) 


Substitution  of  the  strain-displacement  Equation  (2.3)  into  the  constitutive  relation 


given  by  Equation  (2.S)  provides  another  useful  identity  as  shown  by  Equation  (2.1 1). 
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dw  5u 

dx  dz 

(2.11) 


Equation  (2.4)  is  now  substituted  into  the  left  hand  side  column  vector  of  Equation  (2. 1 1). 
Both  Equation  (2.10)  and  the  modified  Equation  (2.1 1)  are  substituted  into  Equation  (2.9) 
to  produce  Equation  (2.12). 
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(2.12) 
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Until  this  point,  the  elasticity  relations  have  been  manipulated  to  establish 
displacement  as  the  variable  for  which  to  solve.  Now,  to  discretize  the  problem  for 
algebraic  computational  solution,  assume  that  over  a  small  domain  each  of  the 
displacements  can  be  represented  by  a  polynomial.  By  discretizing  a  three  dimensional 
domain,  small-but-finite-sized  volume  elements  can  be  formed.  These  finite  elements 
become  eight-noded.  Each  node  can  have  three  orthogonal  displacements.  Then,  each 
element  has  a  total  of  24  degrees  of  freedom  (dof).  Calling  'Uj',  'v,',  and  'w,'  the 
displacement  at  each  node,  displacements  can  be  defined  in  terms  of  polynomial  shape 
functions  as  follows  in  Equation  (2.13). 

“  =EH.u,  ;  V  =5^Hv;  w  =  H,w,  (2  13) 

I  •!  i  1  »1 


Also,  the  first  partial  derivatives  of  the  shape  functions  exist  as  shown  in  Equation 
(2.14). 


i-1  ^ 

53  A  ’ 
1-1  ^ 


i-I 

i-l 

i-l  ^ 


A 

7  - -'^1 

h  4  ' 


(2.14) 


The  '3x1'  displacement  vector  can  be  expressed  as  the  product  of  a  '3  x  24'  shape 
function  matrix  and  a  '24  x  1'  nodal  displacement  vector  as  shown  in  Appendix  A.  The 
product  of  the  '6  x  3'  partial  differential  matrix  of  Equation  (2. 12)  and  the  '3  x  24'  shape 


function  matrix  are  typically  combined  into  one  '6  x  24’  matrix.  This  matrix  is  referred 
to  as  the  'B'  matrix  in  this  development  and  the  full  'B'  matrix  is  also  shown  in  Appendix 
A.  In  shorthand  notation  shown  in  Equation  (2. IS),  the  'B'  matrix  can  be  partitioned  into 
sub-matrices  labelled  'B,',  where  i=l  to  8  for  eight  sub-matrices. 


W  ■[MWP.lWNWNP.il 


(2.15) 


Based  on  this  shorthand  notation,  the  sub-matrices  are  defined  in  Equation  (2.16). 


an 

0  — 10 
dy 


(«.1  = 


dU,  6H 
— 1  — 10 
dy  5x 

dH  aH 

0  — 1  — 1 
dz  dy 

an  aH 

— 10  — 1 
dz  ax 


where  i  =  1  to  8 


(2.16) 


The  Galerkin  method  (or  more  exactly  the  Bubnov-Galerkin  method)  takes  the 


weighting  functions  as  equivalent  to  the  shape  functions  as  shown  in  Equation  (2.17). 
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H, 

H: 

H3 

H. 

H, 

H. 

H3 
H. 

Thus  the  Galerkin  method  only  requires  that  the  weighting  function  be  continuous  over 
small  discrete  intervals  which  correspond  to  the  sides  of  the  finite  elements.  Based  on 
the  Galerkin  weighting  functions  shown  above,  the  old  weighting  function  matrix  of 
Equation  (2. 1 2)  can  now  be  written  in  terms  of  the  shape  functions  as  per  Equation  (2.18). 

aw,  aw,  aw, 

ax  dz 

aw.  aw,  aw, 

0  — ?  0  — I  — I  0 

dy  d\  dz 

aw,  aw,  aw, 

0  0  - 1  0  - 1  - ! 

dz  dy  dx. 

Thus,  Equation  (2.12)  can  be  revised  as  shown  in  Equation  (2.19). 

ilf[BY[D][B]dV{d]=  jj  4>y{H,]\dS  (2.19) 
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In  Equation  (2.19),  {d}  is  the  displacement  vector  which  has  grown  from  the  original  '3 
X  I'  vector  of  'u',  V,  and  'w'  to  a  '24  x  1'  vector  of  terms  'u,',  'Vj',  and  'Wj'.  Also,  the 
right  hand  side  surface  integral  '3  x  1*  vector  has  a  vector  for  each  term  and  is  simply 
shorthand  for  a  '24  x  1'  vector  of  discretized  surface  boundary  tractions.  The  original 
three  partial  differential  equations  of  elasticity  equilibrium  have  now  been  reduced  to  a 
matrix  equation  of  size  24.  The  integrals  in  Equation  (2.19)  may  be  easily  solved  if 
simple  forms  for  shape  functions  are  chosen  and  if  the  modelled  section  geometry  is 
relatively  simple. 

If  the  modelled  section  geometry  is  complex,  the  finite  element  geometry  can  be 
made  simple  by  employing  a  transform  mapping  to  another  reference  space.  In  order  to 
map  'x',  y,  and  'z'  coordinates  for  an  irregular  shaped  element  onto  'r',  's'  and  't' 
coordinates  for  a  rectangular  parallelpiped,  the  Jacobian  transform  matrix  is  required.  The 
Jacobian  transform  matrix  essentially  scales  each  component  of  the  original  space  to  a 
new  space.  Thus,  Equation  (2.20)  shows  the  Jacobian  for  an  'xyz'  system  mapped  to  an 
'rst'  system: 

[Jl- 

In  terms  of  the  shape  functions  and  the  nodal  points,  the  Jacobian  is  expressed  in 
Equation  (2.21)  on  the  following  page. 
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(2.21) 


[J]- 


•  an,  ■  cH,  •  3H, 

E-^-y.  E-^ 


ds 
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•an,  'an,  •  an, 

E  E  a-y.  E  — - 


iW 


dt 


However,  since  the  finite  element  integral  equation  includes  shape  function  partial 
derivatives  with  respect  to  the  'xyz'  coordinate  system,  the  inverse  of  the  Jacobian  is 
required  as  described.  Let  [  J  ]  '=  [  F  ],  where  [  F  ]  is  a  '3  x  3’  matrix.  The  shape 
function  derivatives  with  respect  to  the  'xyz'  system  are  now  expressed  in  Equation  (2.22) 
with  respect  to  the  'rst'  system. 


dH. 


aH. 
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an. 
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an,  an. 
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'  +F 

“as  ”  a 


i  =  1  to  8 


(2.22) 


Equation  (2.22)  is  used  to  calculate  the  [  B  ]  and  [  B  matrices  in  the  'rst'  system. 
Equation  (2.21)  is  used  to  transform  the  volume  differential.  The  finite  element  integral 
is  transformed  below  as  shown  in  Equation  (2.23). 


ffflBYlD]lB]dVtt,y.z){<l)  =  IIJlB]’[D][B]\J\dnrAt){d)  (2.23) 

V  V 


The  transformation  to  the  'rst'  coordinate  system  results  in  simplified  finite  element 
integrals.  The  resultant  transformed  elements  are  termed  isoparametric  elements.  The 


term  isoparametric  refers  to  the  equality  between  the  degree  of  the  equations  for 
transforming  'xyz'  into  'rst'  coordinates  and  the  degree  of  the  shape  functions  for 
estimating  displacement.  The  shape  functions  defined  in  terms  of  'rst*  system  variables 
are  shown  in  Equation  (2.24). 


H,  .1(1  «)(l  «)(!  «) 

H,  .  1(1  «)  (l  -s)  (l  »t) 

H,  =1(1  .r)(l  -s)(l  -t) 

H,  =  l(l  «)  (l  »s)  (l  -t) 

*  /  \  /  \  /  (2.24) 

H,  .1(1  -r)(l  *s)(l  »t) 

H.  ■  l(l  -r)  (l  -s)  (l  *>) 

H,  .1(1  -r)(l  -s)(l  -t) 

H,  -l(l  -r)(l  *s)(l  -t) 

In  order  to  compute  the  integrals  with  a  computer  program,  the  Gauss  quadrature 
of  numerical  integration  is  used.  The  volume  integral  is  redefined  as  the  triple  summation 
from  1  to  the  number  of  integration  points  (NIP)  of  the  integrand  evaluated  at  Gauss 
integration  points  (r^,  Sj  and  t,)  and  multiplied  by  weighting  factors  as  shown  in  Equation 
(2.25). 


NIP  NIP  NIP 

fff[BnD]lB]\J\dlir,s,t){d}‘'£'£j:[B]^[D][B]\J\  W,WjW,  {  d  } 

''y 


(2.25) 


The  results  of  the  numerical  integration  may  vary  over  the  elements  in  the  domain  of  the 
model.  Each  of  these  results  can  be  expressed  as  a  '24  x  24'  matrix  which  is  termed  the 
elemental  stiffness  matrix  [  K,  ]. 
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Recall  the  surface  traction  boundary  conditions  from  the  right-hand  side  of  Equation 
(2.12).  The  integration  of  the  directional  component  of  the  applied  stress  over  the 
element  surface  area  of  the  applied  stress  is  equal  to  the  external  force  applied  to  the 
element.  The  result  of  the  integration  is  a  '24  x  1*  vector  {  F,  }  \^ich  is  equivalent  to 
the  force  in  the  'x',  'y',  and  'z'  direction  at  each  of  the  eight  element  nodes.  Substituting 
these  resultant  terms  into  Equation  (2.24)  completes  the  finite  element  derivation  which 
is  shown  in  Equation  (2.26). 

[K.  ]{d)-(F.  )  (2-26) 

The  finite  element  method  involves  combining  these  element  matrix  equations  with 
given  boundary  conditions  in  order  to  form  one  large  system  of  simultaneous  equations 
for  numerical  solution.  The  computer  programming  methodology  to  achieve  such  a 
solution  is  described  in  the  next  section. 

B.  FORTRAN  COMPUTER  PROGRAM  DE\^ELOPMENT 

Finite  element  numerical  analysis  was  accomplished  using  FORTRAN  computer 
programs.  The  FORTRAN  programs  were  implemented  using  the  modular  programming 
concepts  and  FORTRAN  subroutines  developed  by  Akin  [Ref  21].  The  various 
FORTRAN  programs  used  for  this  study  differed  from  each  other  in  only  two  or  three 
subroutines  out  of  approximately  55  subroutines  to  the  main  program. 

Program  housekeeping  data,  nodal  coordinates,  element-node  connectivity,  boundary 
conditions  and  material  properties  are  read  into  the  program  from  one  input  data  file  via 
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several  standard  Akin  subroutines  (for  example,  INPUT,  INVECT,  INPROP,  and 
APLYBC).  Efficient  storage  of  the  input  data,  interim  calculated  program  data,  and 
calculated  results  is  accomplished  through  the  use  of  semi-dynamic  storage.  Two  large 
data  vectors  are  sized  in  the  main  program  to  store  fixed  point  and  floating  point  data 
respectively.  Matrix  or  array  data  is  converted  for  storage  in  these  column  vectors. 
Computation  of  pointers  to  mark  data  locations  in  the  column  vectors  does  require  some 
program  overhead  in  that  up  to  four  subroutines  (for  example,  PTl,  PT2,  PT3,  and  PT4A) 
are  called  for  this  purpose.  However,  the  efficient  use  of  computer  core  memory  through 
the  semi-dynamic  storage  method  provides  a  significant  advantage  in  the  number  of 
elements  and  nodes  which  can  be  used  in  the  model. 

Subroutine  MODEL  provides  overall  control  for  the  rest  of  the  program.  All 
major  subroutines  are  called  from  MODEL.  Subroutine  ASMMDS  computes  the  element 
stiffness  matrices  using  Gauss  quadrature  integration.  (ASMMDS  is  a  locally  developed 
problem  dependent  subroutine).  Two  integration  points  are  used  in  each  coordinate 
direction  at  the  element  nodes.  This  results  in  eight  independent  integrations  points  per 
element.  ASMMDS  implements  the  subroutines  (that  is,  MATER  and  COMDMX)  which 
implement  the  Kwon  micromechanical  cell  method  to  build  the  composite  smeared 
properties.  Fiber  layup  angle  transformations  are  performed  in  subroutine  STNROT, 
Finally,  ASMMDS  calls  subroutine  STORSQ  to  assemble  the  elemental  stiffness  matrices 
into  the  global  stiffness  matrix. 

After  completing  ASMMDS  computations,  program  control  returns  to  MODEL 
where  boundary  conditions  are  applied,  modifications  are  made  to  the  global  force  column 
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vector  and  the  global  stiffness  matrix  (via  subroutine  APLYBC),  and  the  displacement 
solution  is  obtained  (in  subroutines  FACTOR  and  SOLVE).  FACTOR  uses  the 
Cholesky  method  to  factor  the  system  square  matrix  into  the  product  of  a  lower  triangular 
matrix  and  its  transpose.  SOLVE  use  Cholesky -Gauss  methods  of  forward  and  back 
substitutions  to  complete  the  solution.  APLYBC,  FACTOR  and  SOLVE  are  standard 
Akin  subroutines. 

Post-processing  to  calculate  stress,  strain,  stiffness  reduction  or  damage 
accumulation  is  accomplished  in  problem  dependent  subroutine  MSTRES.  MSTRES  is 
an  original  subroutine  developed  and  modified  as  required  for  this  study.  MSTRES  re¬ 
computes  the  elemental  stiffness  matrices,  and  using  the  displacement  solution  previously 
obtained,  implements  the  Kwon  cell  method  to  compute  the  micro-level  cell  strain  from 
the  composite  macro-level  strain.  The  micro-level  cell  strains  are  used  to  compute  micro¬ 
level  cell  stresses  and  macro-level  composite  stresses.  Stiffness  reduction  and  damage 
criteria  are  applied  based  on  the  micro-level  cell  stresses  calculated  in  this  module. 
Program  output  is  arranged  from  this  module  or  from  other  standard  Akin  models  as 
necessary. 

Program  iteration  to  modify  material  properties  is  required  for  computation  of 
damage  accumulation.  Necessary  iteration  control  is  designed  into  control  module 
MODEL  for  damage  studies  which  are  discussed  later. 
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m.  STIFFNESS  REDUCTION  OF  A  CERAMIC  MATRIX  COMPOSTTE 


A.  OVERVIEW 

In  a  typical  composite,  the  stiffness  or  Young's  Modulus  of  the  matrix  is  much  less 
than  that  of  the  fiber.  In  polymer  matrix  composites  (PMCs),  the  ratio  of  fiber  stiffness 
to  matrix  stiffness  may  be  on  the  order  of  100.  In  ceramic  matrix  composites  (CMCs), 
the  ratio  of  fiber  stiffness  to  matrix  stiffness  is  usually  on  the  order  of  10  or  lower.  The 
ratio  of  fiber  failure  strength  to  matrix  failure  strength  follows  an  analogous  relationship 
as  the  stiffness  rations  for  PMCs  and  CMCs.  Ceramic  matrix  materials  are  typically 
stiffer,  stronger  and  more  brittle  than  polymer  matrix  materials.  However,  PMCs  tend 
to  fail  based  on  the  stress  and  strain  at  the  fiber  level,  while  CMCs  tend  to  fail  based  on 
relatively  lower  stresses  and  strains  at  the  matrix  level  [Ref  22].  Matrix  cracking  (or, 
more  precisely,  matrix  micro-cracking)  is  one  of  the  causes  of  initial  failure  of  all 
composites,  and,  in  particular,  for  CMCs.  When  the  matrix  develops  cracks,  the  adjacent 
fibers  must  carry  additional  load.  Thus,  the  effect  of  matrix  cracking  is  seen  on  the 
macro-level  as  a  reduction  in  stiffness  for  the  composite  structure. 

Much  previous  research  had  been  directed  to  model  the  matrix  cracking 
phenomenon.  Aveston  and  Kelly  [Ref  23]  extended  the  models  of  the  late  I960',  for 
understanding  the  stress-strain  behavior  of  composites  with  matrix  cracks.  They  used 
rigorous  shear-lag  analysis  to  calculate  the  fiber  stresses  when  the  matrix  cracks  in  a 
unidirectional  fiber  composite  plate.  They  implement  a  "strain  energy  versus  crack 
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growth"  approach  to  relate  the  interfacial  fiber-matrix  shear  stress  to  the  unbonded  matrix 
longitudinal  stress  after  matrix  cracking  has  occurred.  Hahn  and  Tsai  [Ref  2S]  develop 
a  bilinear  model  to  explain  stress-strain  behavior  before  and  after  matrix  cracking.  Hahn 
and  Tsai  developed  their  model  with  experimental  results  for  a  [0790°],  glass  fiber,  epoxy 
matrix  composite  loaded  in  a  uni-stress  condition.  Garrett  and  Bailey  [Ref  25]  and 
Parvizi  et  al.  [Ref  26]  tested  cross-ply  glass  fiber,  polyester  matrix  composite  materials 
to  observe  matrix  cracking.  They  validated  some  of  the  earlier  theoretical  work  of 
Aveston  and  Kelly.  They  also  studied  the  effect  of  applied  stress  and  matrix  crack 
spacing.  Reifsnider  and  his  team  [Ref  27]  identified  a  matrix  crack  pattern  in  terms  of 
a  Characteristic  Damage  State  related  to  applied  stress  and  based  on  their  experimental 
work  with  composite  laminates  with  multiple  fiber  layup  angles.  Talreja  [Ref  28,  29] 
uses  a  complex,  rigorous  continuum  mechanics  approach  to  derive  constitutive  relations 
which  essentially  extends  the  work  of  the  authors  cited  above.  However,  Talreja's 
approach  is  based  on  factors  adjusted  by  a  curve-fitting  to  experimental  results.  Talreja's 
method  provides  some  insight  into  the  stiffness  reduction  from  matrix  cracking  but 
implementation  of  his  model  requires  significant  experimental  data. 

This  research  attempts  to  extend  the  prior  studies  by  developing  a  correlation  for 
composite  stiffness  reduction  due  to  matrix  cracking.  The  correlation  is  based  on  results 
from  simple  tensile  testing  of  a  composite  specimen.  The  tensile  testing  provides 
characterization  of  stiffness  reduction  versus  applied  load.  While  this  correlation  approach 
is  not  predictive  from  first  principles  and  basic  material  properties,  it  provides  a 
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methodology  to  predict  global  structural  behavior  of  a  composite  material  to  the  micro¬ 
level  stresses  in  the  fiber  and  matrix. 


B.  EXPERIMENTAL  BASIS 

Wang  [Ref  30]  performed  comprehensive  testing  on  continuous  Nicalon  (SiC)  fiber 
reinforced  calcium  aluminosilicate  (CAS)  composite  systems.  Wang  focused  much 
attention  on  standardization  of  specimen  design  and  testing  methodology  due  to  lack  of 
standards  for  testing  fiber  reinforced  ceramic  composites.  Manufacturer  supplied  data  for 
the  fiber  and  matrix  are  indicated  in  Table  3.1  below. 

TABLE  3.1  CHARACTERISTICS  OF  MATRIX  AND  FIBERS  IN  THE  SiC/CAS 

COMPOSITES 


Properties 

CAS  Matrix 

Nicalon,  SiC 
Fibers 

Elastic  Modulus 
(GPa) 

98 

193.2 

Tensile  Strength 
(MPa) 

124 

2760 

Poisson's  Ratio 

0.255 

Not  Provided 

Wang  determined  the  following  properties  for  a  36%  fiber  volume  SiC/CAS  unidirectional 
composite;  Elastic  Modulus  ranged  from  124  to  131  GPa  and  Poisson's  ratio  ranged  from 
0.29  to  0.30.  These  values  agree  with  results  of  the  micromechanical  model. 

Wang  performed  uniaxial  cyclic  loading  of  unidirectional  [0]  ,  [Oj/OO/Oj],  and 
(OJ/9O3/O3]  composites.  Stiffness  reduction  results  of  these  tests  are  used  for  this  study. 
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C.  MATRIX  STIFFNESS  REDUCTION  CORRELATION 


Kwon  and  Bemer  [Ref.  31]  proposed  the  following  method  to  account  for  the  effect 
of  matrix  cracking  on  the  composite  stiffness.  The  degradation  of  the  matrix  elastic 
modulus  is  assumed  to  be  a  function  of  the  micro-level  matrix  stress  as  shown  in 
Equation  (3.1). 


E“ 


(3.1) 


In  Equation  (3.1),  the  superscript  'm'  represents  a  matrix  property,  the  subscript  'o' 
represents  the  original  undamaged  elastic  modulus,  and  'f  is  a  stiffness  degradation 
function  which  depends  on  the  matrix  ceil  stresses.  A  Weibull  type  distribution  function 
of  the  matrix  cell  von  Mises  stress  values  is  used  to  model  the  stiffness  degradation  of 
the  matrix.  The  stiffness  degradation  Weibull  function  is  given  in  Equation  (3.2). 


f  =( 


exp 


-y 


^vme 


'y» 


if  a"  >  a: 


y» 


(3.2) 


•  r*  TO  ^  TO 


In  Equation  (3.2),  is  the  equivalent  von  Mises  stress  of  the  matrix,  is  the  yield 
or  failure  strength  of  the  matrix,  and  y  and  [3  are  constants  of  the  material  which  are 
determined  from  experiment. 

A  one-eighth  model  of  a  tensile  test  specimen  was  developed  for  input  into  an  FEM 
program  developed  as  described  in  Chapter  II.  The  composite  model  consisted  of  two 
layers  of  20  elements  for  a  total  of  40  elements.  Material  properties  from  Table  3.1,  CAS 
matrix  Poisson's  ratio  of  0.31,  and  fiber  volume  of  36%  were  input  for  the  analysis.  A 
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CAS  matrix  yield  stress  value  of  130  MPa  was  used  for  the  stiffness  reduction  function. 
Boundary  conditions  of  known  displacements  which  were  derived  from  the  Wang 
experimental  strain  and  composite  stiffness  reduction  data  were  also  input  into  the  FEM 
program.  The  FEM  postprocessor  determined  the  matrix  and  fiber  cell  stresses.  The 
matrix  von  Mises  level  stress  values  were  calculated  at  the  cell  level  and  then  averaged 
throughout  the  model.  Average  matrix  cell  level  von  Mises  stress  values  were  considered 
valid  as  the  specimen  was  modelled  for  uniaxial  displacement  loading.  A  separate 
FORTRAN  program  based  on  the  Kwon  micromechanical  model  as  implemented  by 
subroutine  COMDMX  was  run  to  correlate  the  matrix  stiffness  reduction  to  the  composite 
stiffness  reduction.  This  step  essentially  was  directed  at  determining  what  the  required 
level  of  matrix  stiffness  degradation  would  be  for  a  given  composite  stiffness  degradation. 
At  this  point,  necessary  information  to  perform  a  correlation  was  accumulated.  This  data 
is  displayed  in  Table  3.2.  Next,  Equation  (3.2)  was  manipulated  (using  logarithms  and 
basic  algebra)  into  the  format  required  to  perform  a  least  squares  fit  for  constants  y  and 
3-  The  least  square  fit  constants  were  used  to  calculate  the  reduced  matrix  stiffness  in 
order  to  finally  compute  the  new  predicted,  reduced  composite  stiffness.  This  predicted, 
reduced  composite  stiffness  was  correlated  to  the  experimentally  determined,  reduced 
composite  stiffness.  The  least  square  fit  constants  were  adjusted  to  optimize  the 
maximum  error  and  the  root  mean  square  error  between  the  predicted  and  experimental, 
reduced  composite  stiffness  values.  These  steps  were  performed  using  a  Microsoft  Excel 
4.0  spreadsheet,  its  built-in  functions,  and  an  interactive  "macro"  or  program. 
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TABLE  3  2  STIFFNESS  REDUCTION  CURVE  FIT  DATA 


Strain 

%  Stiffness 

E(  Multiplier 

Required  Eg, 

Matrix  von 

Reduction 

Multiplier 

Mises  Stress 

(MPa) 

(experimental) 

(experimental) 

(calculated) 

(computed) 

(computed) 

0.0015 

1.5% 

0.985 

0.968 

147 

0.0020 

4.0% 

0.960 

0.915 

196 

0.0025 

12.5% 

0.875 

0.736 

245 

0.0030 

18.0% 

0.820 

0.620 

294 

0.0035 

22.5% 

0.775 

0.525 

343 

0.0040 

24.5% 

0.755 

0.483 

392 

0.0045 

26.0% 

0.740 

0.452 

441 

0.0050 

28.5% 

0.715 

0.398 

490 

0.0055 

32.0% 

0.680 

0.325 

539 

From  the  above  data  and  analysis,  y  was  determined  to  equal  0.307  and  P  was  determined 
to  equal  1.16.  These  material  constants  were  assumed  to  be  independent  of  the  composite 
layer  configuration,  the  material  properties  of  the  fiber,  and  the  volume  fraction  of  the 
fiber. 


D.  FEM  STIFFNESS  REDUCTION  RESULTS 

The  correlation  for  stiffness  reduction  discussed  above  was  incorporated  into  the 
computer  FEM  program.  To  achieve  stiffness  reduction,  the  program  was  run  once  with 
initial  undamaged  material  properties  to  determine  the  matrix  micromechanical  stress 
values.  The  matrix  stiffness  reduction  factor  was  calculated  based  on  these  stress  values. 
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The  matrix  properties  were  degraded  accordingly,  and  the  FEM  program  was  re-run  to 
determine  the  final  specimen  displacement  and  strain  with  the  matrix  cracking  damage 
present.  Two  different  laminates  were  studied:  Two  cross-ply  composite  laminates  made 
of  SiC/CAS  with  a  fiber  volume  of  36%.  Numerical  results  were  compared  with  Wang's 
experimental  results  as  shown  in  Figures  3.1,  3.2,  and  3.3.  The  numerical  method  for  the 
cross-ply  [O3/90/O]]  composite  predicts  slightly  lower  stiffness  reduction  at  low  composite 
strains  and  slightly  higher  stiffness  reduction  at  medium-to-high  composite  strains  than 
that  obtained  by  experiment.  For  the  cross-ply  [O3/9O3/OJ]  composite,  the  numerical 
method  significantly  under-predicts  the  stiffness  reduction  as  compared  to  the 
experimental  results.  In  the  latter  case,  post-experiment  composite  analysis  indicated 
some  fiber  breakage  in  addition  to  the  matrix  cracking  phenomenon.  The  fiber  breakage 
offers  a  partial  explanation  for  the  higher  than  predicted  stiffness  reduction.  Also,  other 
non-modelled  micro-level  stresses  or  load  transfer  mechanisms  might  account  for  the 
additional  stiffness  reduction. 

To  gain  further  insight  into  the  stiffness  reduction  mechanism  for  the  cross-ply 
composite  laminates,  the  percent  stiffness  reduction  was  plotted  for  matrix  cracking  in 
each  layer  and  for  matrix  cracking  in  the  whole  laminate.  Figures  3.4  and  3.S  display 
these  plots.  One  might  expect  a  greater  degree  of  matrix  cracking  to  occur  in  the  90° 
layer,  as  its  matrix  supports  a  majority  of  the  stress  which  results  from  the  applied 
longitudinal  load  than  does  its  fibers.  But,  for  the  [O3/9O/O3]  composite,  most  stiffness 
reduction  occurred  due  to  matrix  cracks  in  the  0°  layer.  The  reason  is  that  the  0°  layer 
is  much  thicker  than  the  90°  layer.  As  a  result,  the  0°  layer  carries  more  load  compared 
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to  the  90“  layer.  For  the  [O3/9O3/O3]  composite,  the  numerical  method  indicates  that 
stiffness  reduction  due  to  matrix  cracking  in  both  layers  is  comparable.  The  relatively 
thick  90°  layer  in  [O3/9O3/O3]  yielded  more  contribution  to  the  stiffness  degradation,  due 
to  matrix  cracking  in  the  90“  layer,  than  the  relatively  thin  90“  layer  in  [O3/9O/O3]. 

Figures  3.6,  3.7,  and  3.8  display  the  longitudinal  fiber  stress  changes  as  the  matrix 
cracking  is  modelled.  For  all  composite  laminates  considered,  as  the  matrix  cracks,  the 
fiber  stress  is  increased  due  to  load  transfer  from  damaged  matrix  to  the  relatively  stiff 
longitudinal  fibers.  For  the  cross-ply  [03/90/Oj]  laminate  of  Figure  3.7,  matrix  cracking 
in  the  0“  layer  accounts  for  most  of  the  stress  increase  in  the  longitudinal  fibers.  For  the 
cross-ply  [O3/9O3/O3]  laminate  of  Figure  3.8,  matrix  cracking  in  either  layer  contributes 
approximately  the  same  relative  share  of  longitudinal  fiber  stress  increase.  These  fiber 
stress  graphical  results  further  indicate  the  utility  of  the  Kwon  cell  method. 
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Figure  3.1  Experimental  And  Predicted  Stiffness  Reduction  Due  To  Matrix  Cracking  For 
A  Unidirectional  SiC/CAS  Composite  Laminate  Loaded  Uniaxially. 


Figure  3.2  Experimental  And  Predicted  Stiffness  Reduction  Due  To  Matrix  Cracking  For 
A  [03/90/Oj]  SiC/CAS  Composite  Laminate  Loaded  Uniaxially. 
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Figure  3.3  Experimental  And  Predicted  Stiffness  Reduction  Due  To  Matrix  Cracking  For 
A  [O3/9O3/O3]  SiC/CAS  Composite  Laminate  Loaded  Uniaxially. 


COMPOSITE  STRAIN 


Figure  3.4  Stiffness  Reduction  Due  To  Matrix  Cracking  In  Different  Layers  Of  A 
[O3/9O/O3]  SiC/CAS  Composite  Laminate  Loaded  Uniaxially. 
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Figure  3.5  Stiffness  Reduction  Due  To  Matrix  Cracking  In  Different  Layers  Of  A 
[O3/9O3/O3]  SiC/CAS  Composite  Laminate  Loaded  Uniaxially. 


COMPOSITE  STRAIN 


Figure  3.6  Effects  On  Longitudinal  Fiber  Stress  Of  Stiffness  Reduction  Due  To  Matrix 
Cracking  Of  A  Unidirectional  SiC/CAS  Composite  Laminate  Loaded  Uniaxially. 
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Figure  3.7  Effects  On  Longitudinal  Fiber  Stress  Of  Stiffness  Reduction  Due  To  Matrix 
Cracking  In  Different  Layers  Of  A  [03/90/Oj]  SiC/CAS  Composite  Laminate  Loaded 
Uniaxially. 


COMPOSITE  STRAIN 


Figure  3.8  Effects  On  Longitudinal  Fiber  Stress  Of  Stiffness  Reduction  Due  To  Matrix 
Cracking  In  Different  Layers  Of  A  [OJ/9O3/O3]  SiC/CAS  Composite  Laminate  Loaded 
Uniaxially. 
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E.  PARAMETRIC  STUDY 

A  parametric  study  for  different  fiber  volumes  in  a  unidirectional  SiC/CAS 
composite  laminate  was  conducted  and  reported  by  Kwon  and  Bemer  [Ref  32],  Results 
for  stiffness  reduction  sensitivity  to  the  variation  of  fiber  volume  are  indicated  in  Figure 
3.9.  As  the  magnitude  of  the  applied  strain  boundary  condition  was  increased,  the  fiber 
stress  also  increased.  This  resulted  in  an  increased  matrix  stress  too.  When  matrix  stress 
increased,  the  stiffness  reduction  factor  also  increased.  Thus,  composite  stiffness 
reduction  increased  for  all  fiber  volumes  modelled  when  the  strain  was  increased.  At 
higher  fiber  volumes,  the  numerical  model  indicated  smaller  stiffness  reduction.  This 
result  is  expected:  As  the  amount  of  fibers  was  increased,  the  share  of  load  carried  by 
the  fibers  increased  likewise.  Thus,  local  matrix  stress  and  subsequent  matrix  cracking 
was  decreased. 

Figure  3.10  indicates  the  longitudinal  fiber  stress  versus  strain  for  several  different 
fiber  volume  fractions.  The  fiber  stress  increases  with  strain  due  to  increasing  load 
transfer  from  the  damaged  matrix  to  the  undamaged  fibers.  As  the  effect  of  matrix 
cracking  is  reduced  with  higher  fiber  volumes,  the  rate  of  increase  of  fiber  stress  goes 
down. 

The  parametric  study  results  appear  reasonable.  However,  comparison  with 
experimental  data  is  warranted  to  fully  validate  the  results.  Also,  experimental 
verification  could  substantiate  the  assumption  that  the  stiffness  reduction  correlation, 
developed  in  section  'C  above,  is  truly  independent  of  the  fiber  volume. 
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Figure  3.9  Stiffness  Reduction  Versus  Strain  For  Different  Fiber  Volumes  For  A 
Unidirectional  SiC/CAS  Composite  Laminate. 
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Figure  3.10  Longitudinal  Fiber  Stress  Versus  Strain  For  Different  Fiber  Volumes  For  A 
Unidirectional  SiC/CAS  Composite  Laminate. 
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IV.  STIFFNESS  REDUCnON  MODIFIED  BY  THERMAL  RESIDUAL  STRESS 


A.  BACKGROUND 

The  mismatch  between  the  coefficient  of  thermal  expansion  (CTE)  for  the  fiber  and 
the  matrix  of  ceramic  matrix  composites  (CMCs)  and  the  relatively  high  fabrication 
temperatures  for  CMCs  can  cause  significant  residual  thermal  stresses.  These  residual 
thermal  stresses  may  have  a  significant  effect  on  the  performance  of  CMCs. 

The  micromechanical  matrix  stress  was  the  primary  factor  in  the  stiffness  reduction 
correlation  developed  for  the  SiC/CAS  as  discussed  in  the  previous  chapter.  The  Weibull 
type  correlation  has  a  reasonable  physical  basis  in  using  the  ratio  of  the  difference  of 
matrix  cell  stress  and  matrix  yield  stress  to  the  matrix  yield  stress:  This  ratio  indicates 
some  sort  of  threshold  matrix  cell  stress  level  at  which  the  matrix  begins  to  develop 
cracks.  Residual  thermal  stresses  in  the  matrix  could  affect  the  threshold  stress  for  the 
onset  of  cracking. 

A  calculation  was  performed  to  estimate  the  residual  thermal  stresses  in  the 
SiC/CAS  composite  laminates  which  were  analyzed  in  Chapter  III.  The  calculated 
thermal  residual  stress  terms  were  incorporated  into  the  FEM  program  for  the  analysis  of 
stiffness  reduction.  Results  of  the  numerically  predicted  stiffness  reduction  were 
compared  with  experimental  data. 
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B.  DERIVATION 


Two  separate  calculations  were  performed  for  the  thermal  residual  stress  in  the 
matrix  of  the  SiC/CAS  composites:  unidirectional  laminate  calculations  and  cross-ply 
laminate  calculations.  The  coefficient  of  thermal  expansion  (CTE)  for  the  SiC  fiber,  Of, 
was  taken  as  4.0  x  lO"*  per  degree  Celsius;  die  CTE  for  the  CAS  matrix,  a,,  5.0  x  lO"* 
per  degree  Celsius.  This  data  is  based  on  Wang  [Ref.  30].  Other  basic  material 
properties  listed  previously  in  Table  3.1  were  also  used.  A  temperature  delta  of  300‘’C 
was  used  as  a  typical  value  for  the  difference  between  fiber/matrix  bond  setting 
fabrication  temperature  to  application  temperature.  The  approximate  computation 
described  below  relies  on  the  Chamis  relations  [Ref.  3]  to  determine  smeared  properties 
for  the  composite  laminates.  A  method  to  numerically  compute  the  thermal  residual 
stresses  using  the  Kwon  cell  method  is  briefly  outlined  in  Appendix  B. 

1.  Unidirectional  Case 

Consider  a  two  dimensional  view  of  the  fiber  and  matrix  interface  as  two  solid 
blocks  which  are  adjacent  to  one  another.  A  force  balance  must  exist  between  the  bonded 
fiber  and  matrix  after  cooldown  from  bond  setting  during  fabrication  as  shown  in 
Equation  (4.1). 

OfAf=a„A„  (4.1) 

If  each  side  of  Equation  (4.1)  is  multiplied  by  the  composite  thickness  and  then  divided 
by  the  total  composite  volume.  Equation  (4.2)  is  obtained. 

(4.2) 
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The  mechanical  strain  of  the  fiber  and  matrix  is  given  by  the  basic  Hooke's  Law 
relationship  and  is  equal  to  the  respective  elastic  modulus  divided  by  the  respective  stress 
term.  The  thermal  strain  is  equal  to  the  appropriate  CTE  times  the  temperature  delta. 
Equation  (4.3)  shows  that  the  sum  of  the  mechanical  strains  of  the  fiber  and  the  matrix 
must  be  equal  to  the  magnitude  of  the  difference  of  the  thermal  strains  of  the  matrix  and 
the  fiber. 


o,  o 

+_=.  =1  a  -a,  I  AT 
PR  I  ■  M 


Ef 


(4.3) 


When  Equations  (4.2)  and  (4.3)  are  combined  algebraically.  Equation  (4.4)  results. 


,  I  -ctt  I  at 


1.^ 


(4.4) 


Although  developed  independently,  Equation  (4.4)  is  equivalent  to  the  result  derived  by 
Sambell,  et.  al.,  [Ref  33].  When  the  qipropriate  values  are  substituted  into  Equation 
(4.4),  the  resultant  residual  thermal  stress  in  the  matrix  is  approximately  1 S  MPa. 


2.  Cross-Ply  Case 

The  cross-ply  development  for  thermal  residual  stress  is  analogous  to  the 
unidirectional  development,  but  necessarily  more  complex  due  to  the  consideration  of  the 
two  laminae  of  different  fib.r  layup  angle  and  the  two  stress  directions.  In  the  following 
derivation,  subscript  'L'  represents  the  global  longitudinal  direction,  subscript  T'  represents 
the  global  transverse  direction,  superscript  'O'  represents  the  0°  fiber  layup  lamina,  and 
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superscript  '90'  represents  the  90°  fiber  layup  lamina.  Four  different  strain  terms  are 
represented  in  Equation  (4.S). 


Because  the  strains  in  the  global  transverse  and  longitudinal  directions  must  be  equal. 
Equation  (4.S)  actually  represents  two  independent  equations  with  four  unknowns.  To  add 
two  required  equations  to  the  system.  Equation  (4.6)  describes  the  force  balance  between 
the  laminae. 


=0 

A‘’o?M’°or  =0 

To  solve  for  the  longitudinal  and  transverse  CTE  values  and  elastic  modulus  values  for 
the  fiber,  the  equation  of  Chamis  [Ref  3]  were  used.  The  matrix  CTE  was  assumed  to 
be  isotropic.  Noting  that  the  ratio  of  the  areas  for  the  0°  lamina  to  the  90°  lamina  is  two 
for  the  [O3/9OJ/O3]  composite  and  six  for  the  [O3/9O/O3]  composite,  the  system  of  four 
equations  and  four  unknowns  given  by  Equations  (4.5)  and  (4.6)  is  solved  to  give  residual 
thermal  stress  values.  Matrix  thermal  residual  stress  values  are  shown  in  Table  4.1. 
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TABLE  4.1  MATRIX  THERMAL  RESIDUAL  STRESSES 


Lamina 

[O3/9O3/O3]  Laminate 
Matrix  Thermal  Residual 
Stress  (MPa) 

[O3/9O/O}]  Laminate 
Matrix  Thermal  Residual 
Stress  (MPa) 

Oi 

02 

O2 

0“ 

-19.3 

25.3 

-8.43 

10.6 

90“ 

-40.8 

52.1 

-47.3 

68.3 

C  FEM  STIFFNESS  REDUCTION  WITH  THERMAL  RESIDUAL  STRESS 

The  Weibull  stiffness  reduction  function  constants  were  re-computed  using  the 
matrix  residual  stress  values.  The  new  y  was  0.238  and  the  new  P  was  1.42.  Figure  4.1 
shows  the  graph  for  the  stiffness  reduction  of  the  unidirectional  composite  laminate.  The 
new  predicted  stiffness  reduction  values  show  good  agreement  with  the  experimental 
results,  but  the  overall  comparison  indicates  mixed  results  for  this  case.  The  new 
predicted  stiffness  reduction  values  show  slight  improvement  over  the  previous  (without 
thermal  residutd  stress)  predicted  values  at  low  strains.  The  previous  model  shows 
slightly  better  agreement  with  experimental  results  at  medium  to  high  strain  rates.  Both 
ptedictions  are  within  the  same  range  of  accuracy  for  the  constants  and  assumptions 
employed  in  deriving  the  model. 

Figure  4.2  shows  the  graph  of  stiffness  reduction  for  the  cross-ply  [O3/9O/O3] 
composite  laminate.  The  new  predicted  stiffness  reduction  values  show  vastly  improved 
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agreement  with  experimental  results  as  compared  to  the  previous  model.  For  medium  to 
high  strains,  the  new  predicted  values  are  extremely  close  to  the  experimental  values.  For 
low  strain  levels  the  predicted  stiffness  reduction  due  to  matrix  cracking  is  smaller  dian 
both  the  experimental  results  and  the  previous  results.  However,  the  shape  of  the  stiffness 
reduction  curve  shows  upward  curvature  and  then  levels  off  just  short  of  achieving 
downward  curvature.  This  type  of  behavior  might  occur  if  weak  fibers  in  the  composite 
laminate  were  cracked  or  damaged  at  a  low  strain  level.  From  a  statistical  viewpoint, 
some  early  weak  fiber  breakage  is  a  real  probability. 

Figure  4.3  shows  the  graph  of  stiffness  reduction  with  thermal  residual  stress  for  the 
cross-ply  [O3/9O3/O3]  composite.  The  numerical  prediction  with  thermal  residual  stress 
incorporated  shows  slight  improvement  in  matching  the  experimental  results  as  the 
previous  numerical  method.  These  results  indicate  the  need  for  further  refinement  of  the 
stiffness  reduction  correlation  for  thermal  residual  stresses  and  their  effect  on  matrix 
cracking. 
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Figure  4. 1  Experimental  And  Predicted  Stiffness  Reduction  Due  To  Matrix  Cracking  For 
A  Unidirectional  SiC/CAS  Composite  Laminate  Loaded  Uniaxially. 


Figure  4.2  Experimental  And  Predicted  Stiffness  Reduction  Due  To  Matrix  Cracking  For 
A  [0j/90/0j]  SiC/CAS  Composite  Laminate  Loaded  Uniaxially. 
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Figure  4.3  Experimental  And  Predicted  Stiffness  Reduction  Due  To  Matrix  Cracking  For 
A  [O3/9O3/O3]  SiC/CAS  Composite  Laminate  Loaded  Uniaxially. 
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V.  DAMAGE  PROGRESSION  AND  FAILURE  IN  A  POLYMER  MATRIX 


coMPOsrra 


A.  OVERVIEW 

A  finite  element  analysis  was  also  performed  concerning  the  effects  of  flaws  (in  this 
case,  a  circular  hole)  on  the  failure  behavior  of  composite  laminates.  From  classical 
elasticity  based  on  isotropic  materials,  an  "infinite"  plate  containing  a  small  circular  hole 
and  loaded  in  uniaxial  tension,  "o^,"  will  experience  a  stress  concentration  factor  (SCF) 
at  the  edge  of  the  hole.  The  SCF  produces  a  stress  of  magnitude  3  times  in  a  region 
outside  the  hole  and  perpendicular  to  the  applied  stress.  Konisch  and  Whitney  [Ref  34] 
noted  that  the  classical  SCF  approach  produced  anomalous  results  when  applied  to 
composite  laminates  with  holes.  They  modified  a  solution  to  the  elasticity  equations  for 
an  isotropic  material  for  application  to  a  composite  plate  with  a  hole.  Then,  using  a  semi¬ 
inverse  solution  technique,  they  determined  the  stresses  in  the  vicinity  of  the  hole.  Their 
semi-inverse  solution  produced  a  slightly  improved  correlation  with  experimental  data. 
However,  overall  force  equilibrium  was  not  satisfied  throughout  the  composite  by  this 
type  of  semi-inverse  solution. 

The  study  of  stress  around  holes  in  composite  plates  is  predicated  by  the  need  to 
use  bolted  or  pinned  joints  to  connect  such  plates  within  structures.  Bolt  hole  analysis 
has  added  complexity  due  to  the  loading  in  the  hole  area.  Chang,  Scott  and  Springer 
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[Ref.  35,  36]  studied  the  failure  of  such  joints  in  composite  materials.  They  developed 
models  for  the  bolt  hole  boundary  conditions  and  implemented  a  variation  of  the  Yamada 
and  Sun  [Ref  37]  failure  criteria.  The  Yamada-Sun  criteria  is  a  quadratic  equation  for 
composite  laminate  failure  prediction  based  on  average  laminate  properties.  More 
recently,  Chang  and  Chang  have  extended  the  early  work  to  study  the  failure  of  composite 
plates  with  non-loaded  holes  [Ref  38].  The  recent  Chang  research  is  based  upon  FEM 
analysis  using  composite  laminate  theory.  Chang  used  statistical  Weibull  functions  to 
determine  degraded  stiffness  properties  for  the  composite  laminate  when  composite 
damage  criteria  are  exceeded.  Chang  implemented  a  fiber  damage  zone  failure  method 
covered  by  Hahn  and  Tsai  [Ref  39].  The  Chang  numerical  predictions  show  good 
agreement  with  experimental  results.  This  current  work  attempts  to  similarly  predict 
failure  progression  using  the  Kwon  micromechanical  model  with  its  strong  physical  basis 
in  the  micromechanical  properties  of  the  composite  constituents.  Further  finite  element 
analysis  and  modelling  remain  to  fully  compare  this  present  method  with  Chang's  results. 

B.  FAILURE  CRITERIA  AND  MODELLING 

Chang  implemented  failure  criteria  at  the  composite  level  using  modifications  of  the 
quadratic  relations  of  Yamada  and  Sun.  Since  material  failure  has  physical  origins  at  the 
micromechanical  level,  it  is  desirable  to  establish  failure  criteria  on  material  properties  at 
that  level.  Aboudi  [Ref  8]  has  proposed  the  criteria  presented  in  Equations  (5.1)  and 
(5.2). 
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(5.1) 


is.j 


(5.2) 


In  Equation  (5.1),  o',,  is  the  local  fiber  longitudinal  stress  and  X,  is  die  ultimate  tensile 
strength  of  the  fiber  in  the  longitudinal  direction.  In  Equation  (5.2),  0*22  is  the  local 
transverse  stress  of  the  matrix,  a", 2  is  the  local  matrix  shear  stress  in  the  longitudinal- 
transverse  plane,  X„  is  the  ultimate  tensile  strength  of  the  matrix,  and  S.  is  the  ultimate 
shear  strength  of  the  matrix.  Thus,  equations  (5.1)  and  (5.2)  are  based  on  local  cell  stress 
levels  of  the  fiber  and  matrix  and  on  the  ultimate  strengths  of  the  fiber  and  matrix.  Local 
fiber  and  matrix  cracking  damage,  indicating  local  laminate  failure,  is  assumed  to  occur 
when  these  criteria  are  exceeded. 

When  the  composite  level  stresses  exceeded  his  global  failure  criteria,  Chang 
modelled  failure  using  a  different  method  for  both  the  fiber  and  composite.  For  the  fiber, 
Chang  degraded  fiber  properties  by  reducing  the  fiber  transverse  elastic  modulus  to  zero 
and  the  fiber  longitudinal  elastic  modulus  and  shear  modulus  according  to  an 
experimentally  determined  Weibull  distribution.  Chang  degraded  the  matrix  properties 
by  reducing  the  matrix  transverse  elastic  modulus  to  zero.  For  this  work,  when  the  fiber 
or  matrix  criteria  of  equations  (5.1)  and  (5.2)  was  exceeded,  the  respective  elastic  moduli 
and  shear  moduli  for  all  directional  combinations  were  reduced  by  a  factor  to  near  zero. 
Actual  reduction  to  zero  was  not  possible  using  numerical  methods  due  to  resultant 
singularities  in  the  matrix  equations. 
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C  FEM  DAMAGE  PROGRESSION  STUDY 


As  noted  above,  comparison  to  and  extension  of  the  Chang  research  was  a  desired 
goal.  Thus,  the  material  chosen  for  this  numerical  analysis  was  based  on  the  Fiberite 
T300/1034-C  graphite  epoxy  composite  which  was  also  used  for  the  Chang  work  A 
cross-ply  polymer  matrix  composite  was  studied;  a  [0/90],  laminate.  The  fiber  and  matrix 
properties  were  "reverse  engineered"  from  the  smeared  composite  properties  presented  by 
Chang.  The  fiber  and  matrix  properties  were  determined  by  an  interactive  trial  and  error 
method  utilizing  the  Kwon  ceil  model.  Table  S.l  below  shows  the  Chang  smeared 
properties  and  the  calculated  fiber  and  matrix  properties.  Fiber  volume  was  calculated 
to  be  50%. 

TABLE  5.1  GRAPHITE  EPOXY  T300/1034-C  PROPERTIES 


Property 

Composite  - 
Chang  Smeared 

Fiber 

Calculated 

Matrix 

Calculated 

Longitudinal 
Elastic  Modulus 
(Msi) 

21.3 

40.0 

0.50 

Poisson's  Ratio 

0.30 

0.23 

0.40 

Longitudinal 
Tensile  Strength 
(Ksi) 

251 

560 

5.0 

Shear  Strength 
(Ksi) 

19.4 

NA 

23 

Transverse  Tensile 
Strength 

9.65 

NA 

5.0 
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The  FEM  program  discussed  earlier  in  this  report  was  modified  to  simulate  damage 
progression.  The  composite  laminate  was  simulated  with  a  one-eighth  model  due  to 
symmetry.  The  one-eighth  model  laminate  plate  was  0.3  inches  wide,  4.0  inches  long, 
and  0.123  inches  thick.  Hole  radius  was  one-fourth  of  the  laminate  model  width  at  0.123 
inches.  Figure  3.1  displays  a  scaled  version  of  the  specimen  one-eighth  model.  A 
program  control  loop  was  implemented  to  increase  the  applied  end  tension  up  to  72.0  Ksi 
in  increments  of  either  6.0  Ksi  or  7.2  Ksi,  depending  on  the  number  of  load  iterations 
chosen.  The  tension  was  applied  in  the  'x'  direction.  Within  each  load  increment,  a 
conditional  program  loop  tested  for  failure;  recorded  and  wrote  failure  output  data;  and. 
iterated  to  re-compute  degraded  stiffness  properties,  re-solve  for  new  displacements,  and 
re-calculate  to  determine  ceil  stresses  until  the  failure  criteria  indicated  no  further  failures 
for  the  current  load  increment. 

For  the  [0/90],  composite  laminate,  total  matrix  failure  occurred  in  the  90°  layer  at 
approximately  36  Ksi.  The  matrix  failure  occurred  perpendicular  to  the  applied  load 
exactly  along  the  region  of  maximum  stress  according  to  classical  theory.  Total  specimen 
failure  occurred  at  54  Ksi  when  the  fiber  in  the  0°  layer  failed  across  the  specimen.  The 
fiber  failure  around  the  hole  in  the  0°  layer  initiated  78.5°  to  90°  off  the  axis  of  the 
applied  load,  and  then  progressed  along  at  66°  to  78.5°  off  the  axis  of  the  applied  load. 
At  failure,  the  fiber  failed  simultaneously  along  the  outer  specimen  elements  at  66°  to 
78.5°  off  the  axis  of  the  applied  load.  These  results  are  displayed  in  Figures  5.2  through 
5.10  on  the  following  pages.  Chang  modelled  a  [(0/90)5],  laminate  for  numerical 
prediction  of  failure.  For  this  similar  laminate  with  the  same  width-to-diameter  ratio 
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(W/D=4.0),  Chang  observed  failure  at  approximately  52  Ksi  with  the  failure  path  about 
75*  off  the  axis  of  the  applied  load.  Thus,  the  current  results  show  good  agreement  with 
previous  research  in  this  area. 

The  [0/90],  composite  laminate  was  re-modelled  with  end  displacements  applied  as 
compared  to  above,  when  end  forces  which  are  equivalent  to  end  stress  tractions,  were 
applied.  Failure  load  was  within  5%  of  dte  results  above.  The  failure  path  was  also 
similar  to  the  results  discussed  above. 

These  results  differ  only  slightly  from  the  Chang  analytical  results.  A  primary 
cause  for  the  slight  difference  probably  lies  with  the  "reverse  engineered"  fiber  and  matrix 
properties  and  fiber  volume  used  for  this  study.  A  secondary  difference  is  the  failure 
criteria  applied.  The  failure  degradation  model  used  for  this  current  study  was  somewhat 
severe,  in  that  properties  were  reduced  to  the  threshold  of  numerical  singularity. 
However,  the  current  failure  model  is  desirable  from  the  point  of  simplicity. 
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Figure  5.1  -  X-Y  Plane  View  Of  The  FEM  One-Eighth  Model  Of  A  Composite  Plate 
With  A  Circular  Center  Hole. 
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Figure  5.2  -  Damage  Progress  At  12.0  Ksi  Applied  Stress  For  [0/90], 
T300/1034-C  Graphite  Epoxy  Composite  Laminate  Plate  With  Center 
Hole. 


X 


Note:  The  small  circles  indicate  fiber  failure.  The  small  boxes  indicate  matrix  failure. 
Hollow  symbols  indicate  the  0  degree  layer.  Solid  symbols  indicate  the  90  degree 
layer. 
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Figure  5.3  -  Damage  Progress  Al  18.0  Ksi  Applied  Stress  For  [0/90], 
T300/1034-C  Graphite  Epoxy  Composite  Laminate  Plate  With  Center 
Hole. 


Figure  5.4  -  Damage  Progress  At  24.0  Ksi  Applied  Stress  For  [0/90], 
T300/I034-C  Graphite  Epoxy  Composite  Laminate  Plate  With  Center 
Hole. 


57 


X 


Figure  5.5  -  Damage  Progress  At  30.0  Ksi  Applied  Stress  For  [0/90], 
T300/1034-C  Graphite  Epoxy  Composite  Laminate  Plate  With  Center 
Hole. 
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Figure  5.6  -  Damage  Progress  At  36.0  Ksi  Applied  Stress  For  (0/90), 
T300/1034-C  Graphite  Epoxy  Composite  Laminate  Plate  With  Center 
Hole. 
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Figure  S.7  -  Damage  Progress  At  42.0  Ksi  Applied  Stress  For  [0/90], 
T300/1034-C  Graphite  Epoxy  Composite  Laminate  Plate  With  Center 
Hole. 
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Figure  5.8  -  Damage  Progress  At  48.0  Ksi  Applied  Stress  For  [0/90], 
T300/1034-C  Graphite  Epoxy  Composite  Laminate  Plate  With  Center 
Hole. 
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Figure  5.9  -  Damage  Progress  At  54.0  Ksi  Applied  Stress  For  [0/90], 
T300/1034-C  Graphite  Epoxy  Composite  Laminate  Plate  With  Center 
Hole. 
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VL  SUMMARY  AND  CONCLUSIONS 


A.  SUMMARY 

The  new  micromechanical  model  for  composite  properties  developed  by  Kwon  [Ref. 
12]  was  readily  incorporated  into  finite  element  analysis  of  composite  structural 
specimens.  The  Kwon  micromechanical  model  also  provided  for  useful  micro-level  fiber 
and  matrix  stresses  in  the  finite  element  post-processor  subroutines.  The  fiber  and  matrix 
stresses  were  used  to  initiate  stiffness  reduction  correlations  for  material  property 
degradation  and  failure  criteria  for  damage  progression  until  failure. 

The  stiffness  reduction  correlation  developed  for  SiC/CAS  ceramic  matrix 
composites  was  implemented  in  the  finite  element  analysis  based  on  micro-level  matrix 
stresses  determined  from  the  Kwon  micromechanical  cell  model.  A  stiffness  reduction 
correlation  was  developed  by  using  finite  element  analysis  of  a  unidirectional  SiC/CAS 
composite.  The  stiffness  reduction  correlation  was  based  on  a  Weibull  type  function. 
Stiffness  reduction  predictions  for  the  [Oj/90/03]  SiC/CAS  composite  showed  very  good 
agreement  with  experimental  results  of  Wang  [Ref.  30].  Stiffness  reduction  predictions 
for  the  [Oj/OOj/Oj]  SiC/CAS  showed  good  agreement  with  experimental  results. 

Thermal  residual  stresses  were  incorporated  into  the  stiffness  reduction  correlation 
and  finite  element  analysis  for  three  SiC/CAS  composite  laminates.  The  unidirectional 
laminate  was  used  to  determine  the  initial  stiffness  reduction  with  thermal  residual  stress 
correlation.  The  finite  element  analysis  predictions  for  the  [O3/9O/O3]  SiC/CAS  composite 
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showed  markedly  improved  agreement  with  experimental  results.  However,  the  revised 
prediction  for  stiffness  reduction  of  a  [0)/90j/0}]  composite  showed  little  discemable 
improvement  in  matching  experimental  results  as  compared  witfi  the  stiffness  reduction 
correlation  without  thermal  residuals  stress  included. 

Damage  progression  and  failure  studies  were  performed  using  finite  element  analysis 
with  the  Kwon  micromechanical  model  for  a  graphite/epoxy  (T300/1034-C)  cross-ply 
composite  laminate  with  a  center  hole.  The  finite  element  analysis  for  this  study 
predicted  the  failure  load  and  path  with  very  good  correlation  to  the  previous  work  of 
Chang  [Ref  38], 

B.  CONCLUSIONS 

The  micromechanical  cell  model  developed  by  Kwon  [Ref  12]  can  be  readily 
implemented  into  finite  element  analysis  of  composite  structures.  The  initial  studies 
performed  in  this  current  work  indicate  good  correlation  with  experimental  work  and 
previous  research  by  other  authors.  Of  course,  further  analysis  of  other  continuous  fiber 
composite  laminate  structures  with  comparison  to  experimental  results  is  warranted.  Also, 
the  micromechanical  model  need  to  be  modified  for  discontinuous  fiber  composites  of 
whisker  composites. 

The  stiffness  reduction  correlation  developed  in  Chapter  III  needs  further  refinement. 
The  predictions  for  the  [O3/9O3/O3]  cross-ply  laminate  showed  only  fair  correlation  with 
experimental  results.  Incorporation  of  the  characteristic  damage  state  or  zone  concept 
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based  on  matrix  cracking  density  into  the  stiffness  reduction  weibull  function  would 
provide  a  more  clear  relationship  with  the  physical  mechanism  of  the  matrix  cracking. 

Including  thermal  residual  stresses  in  the  stiffness  reduction  correlation  generally 
improved  the  agreement  of  the  predictions  with  experimental  results.  To  further  improve 
these  predictions,  a  full  three-dimensional  thermal  residual  stress  analysis  based  on  the 
Kwon  micromechanical  model  should  be  employed.  Also,  more  detailed  information  on 
composite  fabrication  practice  with  regard  to  heating  and  temperature  might  be  useful  to 
improve  the  computation  of  thermal  residual  stress. 

The  computation  of  fiber  and  matrix  stress  at  the  micromechanical  level  is  an 
exceptionally  powerful  result  of  the  Kwon  micromechanical  cell  method.  This  method 
provides  an  outstanding  physical  foundation  for  future  study  of  composite  structural 
stiffness  degradation,  damage  and  failure  analysis. 
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APPENDIX  A  -  EQUATIONS  FOR  FEM  DERIVATION 
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APPENDIX  B  -  THERMAL  STRESS  DERIVATION  USING  KWON  CELL 

METHOD 


The  following  derivation  has  been  developed  by  Dr.  Young  W.  Kwon.  The 
derivation  is  provided  here  for  completeness  and  reference. 

The  foundation  for  the  Kwon  micromechanical  cell  method  is  based  on  the  relations 
given  in  Equation  (B.l). 

(B.l) 

let 

Based  on  the  unit  cell  as  displayed  in  Figure  1,1,  the  compatibility  equations  for  stress 
and  strain  are  given  in  Equations  (B.2)  and  (B.3). 
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The  total  strain  is  equal  to  the  sum  of  the  mechanical  strain  and  the  thermal  strain 
as  given  by  Equation  (B.4). 
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The  generalized  Hooke's  law  relates  the  mechanical  stress  to  the  mechanical  strain.  From 
Equation  (B.4),  the  mechanical  strain  can  be  composed  of  the  difference  of  the  total  strain 
and  the  thermal  strain.  Equation  (B.5)  shows  this  result. 


An  inverse  relationship  for  the  mechanical  stress  can  also  be  developed  as  shown  in  the 
algebraic  derivation  in  Equation  (B.6)  for  an  orthotropic  material. 
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Equation  (B.7)  re-writes  the  stress  compatibility  equations  in  terms  of  mechanical  strain 
and  thermal  stress  that  were  developed  in  equation  (B.6). 


1. 

oh 

* 

O22 

Eh 

eti 

+E22622  +E23633 

■(”1 

tl2l®ll 

♦E2^E$2  +£^3633 

-« 

2. 

oh 

= 

O22 

E21 

e'l 

■•■E22S22 

+E‘23e‘33 

“2I®II 

+£^2622 

+£23633 

-(“•t 

3. 

oh 

* 

oh 

El, 

s!i 

+E32S22 

+E33E33 

-w, 

I-.C  c 

~  E3,e,| 

■*■£32^22 

+£33633 

-W. 

4. 

= 

at3 

E3. 

et, 

+E^ 

^^32®"  22 

n, 

=  E**  e*' 

^31*^11 

+£32622 

*Eh^h 

-H, 

Equation  (B.8)  re-writes  the  strain  compatibility  equations  in  a  similar  manner. 
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Combining  terms  from  Equations  (B.l),  (B.7),  and  (B.8)  yields  a  system  of  equations 
where  the  mechanical  cell  strains  can  be  expressed  in  terms  of  the  composite  average 
strains  as  shown  in  Equation  (B.9). 
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Equation  (B.9)  can  be  written  in  matrix/vector  equation  format  as  in  Equation  (B.IO). 


Displacement  based  FEM  programs  will  solve  for  composite  macro-level 
displacements.  Postprocessing  modules  first  convert  these  displacement  values  into 
composite  strains.  Inverting  Equation  (B.IO)  allows  solving  for  the  micro-level  cell 
strains.  In  summary,  the  solution  procedure  for  thermal  stresses  becomes; 

1.  Compute  the  smeared  property  matrix  [D]  and  the  thermal  stress  vector  {a,}. 

2.  Compute  the  displacement  vector  {d}  from  [K]{d}={F},  where  vector  (F) 
includes  both  mechanic'll  and  thermal  loads. 

3.  Compute  the  composite  strain  vector{e}. 

4.  Compute  the  fiber  and  matrix  cell  strains{e“}. 

5.  Compute  the  fiber  and  matrix  mechanical  strains  from  the  equation 

+  {e.“}  =  {8“}. 

6.  Compute  the  fiber  and  matrix  cell  stresses  {o“}  from  {o“}  =  [E]  {e„“}. 

7.  Compute  the  composite  stress  {0}  from  {o“}. 
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